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ABSTRACT 

We have developed a time-dependent, multi-energy-group, and multi-angle (S n ) 
Boltzmann transport scheme for radiation hydrodynamics simulations, in one and 
two spatial dimensions. The implicit transport is coupled to both ID (spherically- 
symmetric) and 2D (axially-symmetric) versions of the explicit Newtonian hydrody- 
namics code VULCAN. The 2D variant, VULCAN/2D, can be operated in general 
structured or unstructured grids and though the code can address many problems in 
astrophysics it was constructed specifically to study the core-collapse supernova prob- 
lem. Furthermore, VULCAN/2D can simulate the radiation/hydrodynamic evolution 
of differentially rotating bodies. We summarize the equations solved and methods incor- 
porated into the algorithm and present results of a time-dependent 2D test calculation. 
A more complete description of the algorithm is postponed to another paper. We high- 
light a 2D test run that follows for 22 milliseconds the immediate post-bounce evolution 
of a collapsed core. We present the relationship between the anisotropies of the over- 
turning matter field and the distribution of the corresponding flux vectors, as a function 
of energy group. This is the first 2D multi-group, multi-angle, time-dependent radia- 
tion/hydro calculation ever performed in core collapse studies. Though the transport 
module of the code is not gray and does not use flux limiters (however, there is a flux- 
limited variant of VULCAN/2D), it still does not include energy redistribution and 
most velocity-dependent terms. 

Subject headings: supernovae, multi-dimensional radiation hydrodynamics, transport, 
neutrinos 



1 Racah Institute of Physics, The Hebrew University, Jerusalem, Israel; eli@frodo.fiz.huji.ac.il, ita- 
mar@saba.fiz.huji.ac.il 

2 Department of Astronomy and Steward Observatory, The University of Arizona, Tucson, AZ 85721; bur- 
rows@as.arizona.edu, rwalder@as.arizona.edu 

3 Astronomy Department and Theoretical Astrophysics Center, The University of California, Berkeley, CA 94720; 
thomp@astro.berkeley.edu 

4 Hubble Fellow 



- 2 - 



1. Introduction 

It is now clear that if stellar core collapse simulations are constrained to be spherically- 
symmetric the cores do not supernova (Burrows, Hayes, and Fryxell 1995; Mezzacappa et al. 2001; 
Rampp and Janka 2000,2002; Liebendorfer et al. 2001,2003; Buras et al. 2003; Thompson, Bur- 
rows, and Pinto 2003). However, some calculations performed in 2D (Herant et al. 1994; Burrows, 
Hayes, and Fryxell 1995; Fryer et al. 1999) and in 3D (Fryer and Warren 2002,2003) do explode, 
while others performed in 2D do not (Mezzacappa et al. 1998; Buras et al. 2003). Even though a 
consensus is emerging that multi-D effects are crucial to the mechanism of core-collapse supernovae, 
no censensus has yet emerged concerning which multi-D effects are crucial. The enhancement of 
the efficiency of neutrino heating behind the stalled shock wave due to neutrino-driven convec- 
tion seems to play a role, but this conclusion too remains controversial. It is suspected that the 
approximations to neutrino transport incorporated into the codes which see explosions may have 
compromised the validity of those results. Furthermore, the degree of neutrino heating in the 
semi-transparent region behind the stalled shock may depend upon a proper treatment of multi-D 
transfer in the unstable and convective core. Moreover, the viability of doubly-diffusive instabilites 
in protoneutron stars hinges upon neutrino transport in non-radial directions (Mayle and Wilson 
1988; Bruenn and Dineva 1996). Hence, there is now an emphasis in supernova theory circles on 
improving the treatment of neutrino radiation transport in two and three dimensions. We have 
embarked upon such a program of improvement and this paper is the first in a series on our evolving 
efforts. 

In many astrophysical problems matter and radiation are strongly coupled and are described by 
the system of non-linear radiation hydrodynamics equations. When the optical depth is large every- 
where, the diffusion approximation can be used. In the diffusion limit, the radiation field is in local 
thermodynamic equilibrium (LTE) with the matter and the radiation is described by the non-linear 
diffusion equation. However, when important subregions of the system are optically thin, a full 
Boltzmann problem must be solved. A typical example of the latter is the core-collapse supernova 
problem, where the radiation field consists of several neutrino types, the optical depth varies rapidly 
with time and position, and the coupling to matter may determine the failure or success of the 
neutrino-driven explosion mechanism. In that problem, the radiation field depends upon position, 
direction, frequency (energy), and species, for which any discretization scheme becomes very large 
and expensive. For this reason, to date most time-dependent numerical radiation/hydrodynamics 
conducted in the astrophysical context have been limited to one spatial dimension, or if multi- 
dimensional, has been either flux-limited and gray ("one-group") (Stone, Mihalas, Norman 1992; 
Hayes and Norman 2003) or carried out along radial rays (Rampp and Janka 2000,2002; Buras et 
al. 2003), or both (Burrows, Hayes, and Fryxell 1995). The 2D calculations of Burrows, Hayes, 
and Fryxell (1995) incorporated flux-limited, gray transport along radial rays, and ignored trans- 
verse transport. The 2D calculations of Buras et al. (2003) used a Boltzmann solver, but were 
conducted in ID along various radial rays. Their hydrodynamics is, however, done in 2D. These 
authors correct for transverse radiation pressure and transverse lepton transport using a further 
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update split off from the radial transport. The error introduced by such a procedure is unknown. 
The 2D calculations of Mezzacappa et al. (1998) actually employed the Multi-Group-Flux-Limited 
Diffusion (MGFLD) neutrino heating rates from a previous ID simulation. These deposition rates 
were incorporated into a 2D hydro simulation, but without the feedback that would have arisen 
from the large differences between the ID and 2D density and composition profiles and the shock 
history, as well as from the angular variation seen in 2D. The ZEUS-2D radhydro code of Stone, 
Mihalas, and Norman (1992), and its update by Hayes and Norman (2003), solve the zeroth- and 
first-moment equations of the gray transport equation and use a short-characteristics solution of the 
static transfer equation to obtain the second moment for closure. These codes are implementations 
of the variable Eddington factor technique and avoid some of the pitfalls of flux limiters, but since 
they are gray and not multi- group they are of limited utility for modern supernova simulations. 
The recent increases in computing power now make possible the development of more sophisticated 
codes to perform 2D, multi-angle, multi-group, multi-species radiation hydrodynamics simulations. 

The vast nuclear engineering literature contains numerous methods for solving the Boltzmann 
equation (Adams and Larsen 2002). However, these may not be appropriate for astrophysical 
problems. In astrophysics, we confront stellar configurations in which the optical depth varies by 
many orders of magnitude from their centers to their peripheries and the density and mean free 
paths change rapidly with radius. In this paper, we present a test core-collapse simulation using 
a new algorithm to solve the multi-group Boltzmann equation in one and two spatial dimensions 
which we incorporated into the broader framework of radiation hydrodynamics. We defer a more 
complete discussion of the solver, stencils, and difference scheme to a later work (Livne et al. 2004). 

We summarize aspects of our Boltzmann transport scheme in one and two spatial dimensions 
in §2 and §3. For a full radiation/hydrodynamics code, one needs to couple the radiation field 
to matter through emission, absorption, scattering, and radiation pressure. The resulting time-, 
space-, energy- group-, and angle-dependent code has both ID and 2D (axially-symmetric) variants. 
Section 3 is devoted to a precis of the radiation hydrodynamics scheme in which the new transport 
module is incorporated. We also touch on various subtleties concerning 2D gravity and momentum 
conservation and the importance of freeing the inner core to move along the symmetry axis. In 
§4, we summarize various computational issues concerning problem size and parallelism philoso- 
phy. In §5, we first discuss the results of our comparison of ID stellar collapse simulations using 
VULCAN/ID with the Feautrier /tangent-ray code "SESAME" 5 (Burrows et al. 2000; Thompson, 
Burrows, and Pinto 2003) and then go on to describe our preliminary 2D simulation results in the 
core-collapse context. This test simulation was carried out for 22 milliseconds of physical time. 
We focus on the relationship between an anisotropic matter distribution and its corresponding 
anisotropic radiation field. This is the first 2D time-dependent, multi-group, multi-angle radia- 
tion/hydro calculation ever performed in core collapse studies. For these first test simulations, only 
electron neutrinos are incorporated, though the code is set up to handle the other species. The 2D 
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code can employ unstructured fixed or moving meshes and incorporates the centripetal terms due 
to rotation about the axis of symmetry. Finally, in §6 we list the accomplishments represented by 
this paper, as well as the significant work that remains yet to be done. 



2. The Transport Equations and Approach 

The transport equation for the specific intensity (I) 

-^r + Q-VI + al = S (1) 

c Dt 

describes for all directions fi the additions to and subtractions from a beam of massless particles 
due to the physical processes of absorption, scattering, emission, and streaming. In the notation of 
Morel et al. (1996), I(r, il,e u ,t) is the angular specific intensity (in energy units), a a (r,e u ) is the 
absorption cross section times the number density (inverse absorption mean- free-path), <r s (r, e„) is 
the inverse scattering mean-free-path, a equals a a + a s , S (the source term) equals S em (r, e v ) + a s J 
(where S em is the emissivity), and J(r,e u ,t) equals 4- J Idfl (the zeroth moment of I(r, fl, £ u , t)). 
Scattering is currently assumed to be isotropic, though this is not an inherent feature of our 
formulation, and we use the transport cross section, (1 — {cos9))a^, instead of the total scattering 
cross section (cr^). This approach has a long pedigree and has been shown to work quite well in 
the ID core-collapse context (Burrows et al. 2000; Thompson, Burrows, and Pinto 2003). In eq. 
(1), an energy-group subscript for each energy e u is implied. We solve eq. (1) implicitly and for 
every group and angle. Note that the term \ is the total, Lagrangian time derivative. In eq. (1), 
the other velocity-dependent terms, such as the Doppler (radiation compression) and aberration 
terms (Mihalas and Mihalas 1984), have been dropped. Such velocity effects are not dominant in 
the core-collapse context, but the compression term in particular would need to be included in a 
future generation of supernova simulations if neutrino viscosity effects are to be properly treated. 
Coupling as they do energy groups, angular bins, and spatial zones in the type of global solution 
required for a fully implicit treatment, these two velocity-dependent terms represent one of the 
most complicated numerical challenges of multi-D transport. One component of this challenge, 
as yet unsolved, is to incorporate non-monotonic velocity fields, with the resulting complicated 
characteristic structures, into an implicit solver. This fact is why to date no time-dependent multi- 
group, multi-angle astrophysics code has properly included them. We plan to do so in a later version 
using an operator-split, but explicit, approach that we have developed for the spherically-symmetric 
code SESAME (Thompson, Burrows, and Pinto 2003). 

In spherical symmetry, the one-dimensional transport equation (eq. 1) is 

1 Dl dl (1 - fi 2 ) dl T n 
c Dt or r da 



The discretization of the radial transport equation is easiest when eq. (2) is written in conservative 
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form: 

I PI 1 g(rV) 1 3((1-// 2 )J) _ r 

_ TTT + ~2 — a 1 a h cr7 = 6 . (3) 

c Dt r z Or r Ofi 

This form expresses the fact that the transport operator ft ■ V is the divergence of the flow in 
direction ft. Straightforward integration of this differential equation over a (small) finite volume 
(27rr 2 drdfi) provides a conservative scheme in which appears creation and extinction terms, along 
with the boundary fluxes. 

In 2D axisymmetric geometry, the transport equation we solve becomes 

1 Dl dl 81 81 T n 

cDt + »d-r +r] d-z-^ +C Tl = S (4) 

or in conservative form : 

_ 7T7 H o tt— + aI = S. (5) 

c Dt r or Oz r o<p 

Here, the directional cosines are just the projections of the unit vector ft (in the neutrino/photon 
velocity direction) on the Cartesian axes (Figure 1): 



rj = cos(0) 
\i = sin(0) cos(</>) 
^ = sin(6>) sin(^) , 

where 6 is the polar angle and (ft is the azimuthal angle. Generally, we distribute rj evenly from 
-1 to 1 and (f) evenly from to 7r, but make the number of (j)s a function of r\ in order to tile the 
hemisphere more or less uniformly. We need to calculate only on a hemisphere, and not the entire 
sphere, because of the imposition of axial symmetry. This is also the reason is bounded by and 
7r, and by not 2tt. 

Our numerical scheme uses two features of the transport operator ft ■ V. First, being a first- 
order differential operator along the direction ft, the transport equation can be integrated along 
a characteristic that traverses the computational domain. Second, this operator is the divergence 
of the velocity w v = eft, so that eq. (1) can be regarded as a continuity equation that reflects 
a conservation law for the number of particles. Thus, our aim was to construct a conservative 
difference scheme, unlike the method of short characteristics (e.g., van Noort, Hubeny, and Lanz 
2002), that is motivated in part by the prevalence in astrophysics of radiation hydrodynamics 
problems for which the coupling between the radiation field and matter is stiff. In these cases, 
small relative errors in the energy budget can cause big errors in the dynamics. The core-collapse 
supernova problem, in which the neutrino radiation field contains a non-trivial fraction of the energy 
of the system, is a good example of this type of physical situation. 

In VULCAN/2D, we discretize the angular variables using the discrete-ordinates (S n ) method, 
where the specific intensity depends on a set of representive directions (points on the unit sphere) 
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(Mezzacappa and Bruenn 1993). Multi -group S n methods can be constructed in an efficient way 
for general grids in one, two, or three spatial dimensions. The salient characteristic of S n methods 
is the representation of the angular variable by a set of direction cosines on the unit sphere. Unlike 
other methods, such as the tangent-ray method, the parametrization of the angular variable is 
independent of position and the same set of directions applies for the entire domain. This feature 
makes it easy to generalize the S n method to any dimension and to any structured or unstructured 
grid. A drawback of the S n method is its tendency at low optical depths to concentrate flux along 
the rays. However, in the core-collapse context it is the radiation field interior to the shock wave 
that is most germane, and this region is at low to modest optical depths where the anisotropy of 
the radiation field is not severe. In this work, we limit the discussion to the ID and 2D cases only, 
since the computing resources needed for an accurate 3D solution to the time-dependent transport 
equation are currently unavailable. 



3. The Radiation Hydrodynamics Scheme: Matter/Radiation Coupling, 
Parallelization, and an Aside on Gravity 

The numerical scheme used in VULCAN consists of a Lagrangian step, followed by a remapping 
step to the Eulerian grid (Livne 1993). This makes the code similar in this regard to traditional 
ALE (Arbitrary-Eulerian-Lagrangian) codes. The hydrodynamical variables are all cell-centered, 
except for the position and the velocity, which are node-centered. The variables of the radiation 
field are also cell-centered, so that the interaction between the radiation field and matter is properly 
centered. For radiation hydrodynamics, it is convenient to define the specific intensity as the energy 
specific intensity, rather than the particle specific intensity. Thus, I g is the energy intensity in group 
g (averaged over the group) and I g /su is the particle intensity in group g (averaged over the group), 
e v being the average particle (e.g., neutrino) energy of group g. 

The time advancement in both the radiative and the hydrodynamical sectors is computed in 
the Lagrangian step, whereas the remapping step changes only the spatial discretization of the 
variables over the numerical grid. Therefore, we describe here only the Lagrangian scheme. The 
transport equation itself, for a given source, is always computed in a fully implicit manner. 

We first advance the velocity by half a timestep: 



v n+i/2 = V n + Q 5A ^_W! _ VU G ), (6) 



p n 



where Ug is the gravitational potential (but see the end of this section). The position vector is 
then advanced using 

r n+l =r n + Atv n+l/2 _ ^ 

Denoting by V the volume of a cell, Lagrangian mass conservation takes the form: 

yn 

P U+1 = • (8) 
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We then solve the adiabatic energy equation for the specific internal energy: 

e * = e "-l/2(/+p")(-^-l). (9) 

Equation (9) is iterated to convergence. At this stage, we compute new cross sections and emission 
sources 

S ern = J2 a 9 J 9 q > ( 10 ) 

9 

where Jg 9 is the LTE intensity (a function of density, temperature, and composition). Using 
those cross sections and sources, we solve the transport equation for according to the scheme 
described in Livne et al. (2004). The net change in the radiation energy density is given by: 

AE r = AtJ2^ a g (J e g q -J^ +1 ) (11) 

9 

and this is also minus the net change in the matter energy density. Consequently, we can compute 
the final energy density, pressure, and temperature using 

g n+l = e *_ AEr / p n+l ( 12 ) 

and the equation of state. 

For the supernova problem, we also need to compute the degree of neutronization of matter 
due to electron capture and other charged-current processes (Burrows and Thompson 2003). We 
obtain 

yn + l = yn _ At[ J2 a a {J e q _ jn+ly^ > (13) 

9 ^ aP 

where N a is Avogadro's number and Y e is the elctron fraction. Finally, we advance the velocity due 
to the new matter pressure by a further half timestep and due to the radiation pressure (F"°^ e ) by 
a full timestep: 

v n+i = v n+i/2 + o.5At[JL(-Vp n+1 + 2F™f ) - VU G ] ■ (14) 

The radiation force at grid nodes is evaluated by a simple averaging process using the radiation 
force at cells and the definition of the radiation flux: 

F ™i = \ E / ^ndn . (15) 

C 9 J 

Since we are focussed here on doing radiation/hydrodynamics in two spatial dimensions, not 
three, for reasonable sized grids and EOS and opacity tables we can fit an entire hydrodynamic 
calculation on a single processor. Because the transport solve consumes far more CPU than the 
hydro updating, this suggests a path to parallelism that is both simple to implement and almost 
perfectly scalable. Given that we will handle energy redistribution explicitly, the energy groups are 
not coupled during a timestep. Hence, we can parallelize in energy groups (and neutrino species) and 



-8- 



calculate the hydro on ALL processors concurrently. Each processor handles one energy group for a 
given neutrino species and redundantly calculates the hydro. Interprocessor communication merely 
collects the transport results for the neutrino energy deposition and radiation pressure calculations 
that require integrals over the energy groups and angles. The result is then redistributed over 
the processors at the next timestep and the process is continued. No domain decomposition is 
required. For example, if we follow 20 energy groups and three neutrino species, only 60 processors 
are required. We use MPI to handle the parallelization. In practice, the communication overhead 
is only 2% to 8% of the total run time. 

Gravity is a key force in multi-dimensional astrophysical hydrodynamics. However, many 
calculations in the past have employed only the monopole term and/or have complemented the 
gravitational force term (Fq), written as a gradient of a potential in the momentum equation 
(eq. 6) with a corresponding v ■ Fq term in the energy equation. The latter approach is perfectly 
reasonable, but, given the inherently approximate nature of finite-difference realizations of the 
partial differential equations, does not guarantee momentum conservation nor consistency between 
the momentum and energy equations when written in Eulerian form. To address this, we have 
implemented a version of the code in which the z-component of gravity appears as the divergence 
of a stress tensor (Shu 1992; Xulu 2003). This ensures, in principle, the conservation of momentum 
in that direction, or at least guarantees that in fact the gravitational force of mass parcel "A" 
on "B" is equal and opposite to the gravitational force of mass parcel "B" on "A." This has not 
been the case in past grid-based astrophysical codes that have addressed the supernova problem. 
This approach works with any algorithm employed to calculate the gravitational potential (we use 
a multipole expansion with ~20 terms) and it, along with a gridding scheme such as is described 
in Ott et al. (2004), allows the inner core to move on the grid. All previously published grid- 
based results in supernova theory treated the inner core either in ID (radially) or as a boundary 
of finite radius. Surprisingly, no previous such publication of grid-based results allowed the core 
to move, something that is necessary if we are to address the pulsar kick problem and to treat 
multi-dimensional hydrodynamic effects with realistic gravitational feedbacks. 



4. Computational Issues of Multi-dimensional Problems 

The central computational difficulty in solving the equations of radiation/hydrodynamics in 
many dimensions stems from the fact that, at the very best, the number of operations and the 
number of variables scales with the product of the number of generalized zones in the various 
dimensions involved. VULCAN/2D treats two spatial dimensions (axially-symmetric) , one time 
dimension, one energy-group dimension, and two angular dimensions {9 and <j>). This amounts to a 
total of 6 dimensions. In fact, the number of angular dimensions is the same as would be required 
in three spatial dimensions, so a full 3D simulation would involve only one more dimension and one 
more hemisphere in momentum space (§2). 

In VULCAN/2D, the number of spatial cells (N e ) for a low-resolution simulation might be 
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100 x 100 ~ 10, 000 and for a higher-resolution simulation might be 300 x 300 ~ 100, 000. For an 
S n specific intensity calculation, the number of 9 directions might be n = 10 — 20 and the number 
of 4> directions might average 10, for a total of ~ 15 x 10 = 150 angular bins. The minimum 
number of energy groups required for a good core-collapse simulation is N g = 10 — 20 for each 
neutrino species (Thompson, Burrows, and Pinto 2003). This results in N e x N g x = 15,000,000 
to 300,000,000 specific intensity values at each timestep and for each neutrino type. The radiation 
solve is always much more costly than the hydro solve and completely dominates any run. Doing 
the radiation solve for only an inner fraction of the spatial grid has its merits and can increase 
speed, but we have not yet implemented this. Note that the cost of a corresponding 3D simulation 
would be larger by ~100 zonesx2 hemispheres. 

Currently, on a standard fast processor and for the time-dependent test run described in §5.2 
we required ~0.5 - 2.0 minutes per energy group per timestep. The time-dependent run employed 
the equivalent of 16 x 10 angular bins and 102 x 51 spatial zones. Since we parallelize in energy 
groups and neutrino species (§3) and can expect very good scalability on most clusters, this is a 
benchmark number for the current scheme. Accuracy and stability requirements currently impose 
a timestep after core bounce of ~0.1-0.3 microseconds. Decreasing the number of angular bins (A^) 
can increase the speed dramatically, but with a corresponding loss of accuracy. It is the scattering 
angular redistribution term in eq. (1) that poses the greatest numerical challenges. A problem that 
involves only absorption and emission, or that does not involve large scattering optical depths, is 
very much simpler and faster. 



5. VULCAN Simulations 

5.1. ID Comparison Test: VULCAN/ID versus SESAME 

Though not the focus of this paper, but as a run up to the VULCAN/2D runs, we compared 
collapses using VULCAN/ ID and the ID Lagrangian radiation/hydro code SESAME, based on 
Feautrier variables and the tangent-ray method, as described in Burrows et al. (2000) and Thomp- 
son, Burrows, and Pinto (2003). The progenitor model used for this comparison was the core of 
an 11 M model from Woosley and Weaver (1995). The cross sections were obtained by interpo- 
lating in a three-dimensional table in T, p, and Y e , generated using the formulations of Burrows 
and Thompson (2003). The effects of Pauli blocking by final-state electrons and of stimulated 
absorption, where appropriate, are included. However, only electron-type neutrinos and processes 
were incorporated for these tests. The equation of state we use is the finite-temperature liquid-drop 
model of Lattimer and Swesty (1991) which contains nucleons, NSE nuclei, alpha particles, photons, 
and electrons/positrons, the latter corrected for arbitrary degeneracy and relativity. The spherical 
(ID) calculations were done with the same number of radial zones (300) and energy groups (20). 
The latter span an energy range up to 320 MeV. In the VULCAN/ID calculations, 15 S n angles 
were used, equally-spaced in cos(#). 
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Small differences in zoning and in time stepping cause slight time shifts during the collapse 
phase and there is a 20-millisecond shift between the two calculations in the time of bounce 
(SESAME: 213 msec; VULCAN: 195 msec). The entropies at bounce in the inner 0.5 M dif- 
fer by 0.1 to 0.2 units. The Y e profiles show a corresponding difference, with the values of Y e 
differing by as much as 0.01 units (inside) to 0.05 units (at ~0.6 M Q ) in the inner 1.1 M Q . The 
outgoing deleptonization wave and the Y e trough look quite similar and the overall lepton losses to 
infinity in the two calculations are within 5%. 

From the center to an interior mass of 1.0 M Q , the luminosities for the VULC AN/ID and 
SESAME calculations were within ~1%, but exterior to 1.1 M they differed by up to ~10%. 
This slight discrepancy in the outer regions can be traced to the inherent differences/strengths of 
Eulerian and Lagrangian codes, the different stencils, the inherent differences between tangent-ray 
and S n methods, and differences in the distribution of radial zones after bounce. Interior to ~1.25 
M Q , the Eddington factors for the two runs were the same to within better than a percent, but since 
VULCAN/ ID uses the S n method and S n can't resolve a strongly forward-peaked specific intensity 
distribution without many, many more angular bins, in the outer transparent zones the Eddington 
factor in the VULCAN/ID run saturated at ~0.9. The Eddington factors for the SESAME run 
achieved the correct value of 1.0. However, given that the neutrino-matter energy deposition in the 
so-called "gain" region occurs only in the inner semi-transparent regions, VULCAN/ID, despite 
its S n method, provides in that critical region quite reasonable results. Nevertheless, updates to 
VULCAN that are still necessary are the energy redistribution due to inelastic neutrino-electron 
scattering and the Doppler compression term. 

The differences between the quantities in the ID core-collapse simulations recently compared 
in Liebendorfer et al. (2003) provide a context for our ID comparisons here. The ORNL(Agilc- 
Boltztran)/Garching(Vertex) codes correspond well qualitatively, and in the main quantitatively. 
No evolutionary phases are importantly discrepant. However, at a given time entropies differ by a 
few xO.l, flux factors by 0.05 to 0.1, luminosities by 10-20%, shock radii by 10's of kilometers, and 
velocities by 10's of percent. In addition, there seem to be small numerical issues with all codes 
arising from energy conservation, shock capturing, and the handling of neutrino-matter coupling, to 
name but a few. Nevertheless, despite the different approaches, be they represented by SESAME, 
VULCAN/ID, Agile-Boltztran, or Vertex, ID core-collapse simulations give very similar evolutions 
and results. Note that all these codes indicate that the naively ID (spherical) models do not 
explode. 

5.2. Two-dimensional, Multi-Group, Multi- Angle Simulation Results with 

VULCAN/2D 

There is little in the supernova literature that addresses the issue of aspherical neutrino radi- 
ation fields in 2D (or 3D) matter distributions. Most of it is concerned with the effect of rotation 
on the flux contrast between the pole and the equator and none of it involves the solution of the 
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transport or Boltzmann equation. None of it is multi-group or multi-angle. LeBlanc and Wilson 
(1970) performed 2D flux-limited MHD calculations to investigate the possible role of magnetically- 
driven axial jets during core collapse, but did not address the question of the associated anisotropic 
radiation field. They found that with very rapid rotation an under-energetic MHD jet punched 
along the poles. Symbalisty (1984) too was concerned with the role of magnetic fields during 2D 
collapse, but handled the neutrinos with a crude leakage scheme. Janka and Monchmeyer (1989a,b) 
studied the possible ambiguity in the derivation of the total neutrino burst energy of SN1987a due 
to the anisotropy in the neutrino flux distribution that would be caused by rapid rotation of the 
protoneutron star left behind. They derived a 6 dependence to the neutrino flux from an oblate 
spheroid by associating it predominating with the solid angle subtended by the spheroid at the 
observer. Pole-to-equator contrasts of as much as three were obtained, but no multi-D transport 
calculations was performed. 

Recently, Kotake, Yamada, and Sato (2003) and Shimizu et al. (2001) used a similar approxi- 
mate approach to calculate the flux asphericity for 2D pure- hydro rotating collapse models. They 
used the code ZEUS-2D, but without including the gray, flux-limited transport it contains. They 
estimated the position of the neutrinosphere, derived a neutrino flux enhancement at the pole, and 
concluded that such an asymmetry may play a role in the neutrino-driven mechanism of supernova 
explosions and in the distribution of the debris. We think that there may be merit in this possibil- 
ity and are planning to calculate the anisotropy of the neutrino fluxes for rotating cores using the 
machinery of VULCAN/2D 6 . 

The calculations of Fryer and Heger (2000, 2D) and Fryer and Warren (2003, 3D) are the most 
complete explorations to date of the effects of rotation on the supernova mechanism. However, 
these authors used a flux-limited, gray diffusion approach and did not publish anything concerning 
the derived anisotropies of the radiation field, nor on the potential role of such anisotropies in the 
explosion. Their focus was mainly on the evolution of the angular momentum and on its interaction 
with the convective motions driven by neutrino heating. Finally, as mentioned in §1, Buras et al. 
(2003) conducted state-of-the-art 2D supernova simulations, but used multiple ID radial Boltzmann 
calculations in lieu of real 2D transport. Our VULCAN/2D developments are being undertaken 
in part to improve upon such calculations by including the second angular dimension into full 
multi-group simulations performed in the 2D radiation/hydrodynamic context. 

Our first 2D time-dependent radiation/hydro runs with VULCAN/2D, those that we publish 
here, are to explore aspects of the character of 2D transport in the core collapse context. In a 
sense, they are intended as a "shakedown cruise" of the code. They are not meant to cover the 
full evolution from collapse, through bounce, shock formation, and explosion (?). Here, we study 
the asymmetries of the radiation field using the multi-group, multi-angle code and explore the rela- 



6 Rotation itself, and its effect on the hydro-dynamic flow along the poles (not just the consequent enhancement of 
the driving neutrino fluxes), is also thought to be a potentially important feature in the mechanism of core-collapse 
supernova explosions (Burrows, Ott, and Meakin 2003). 
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tionship between 2D matter distributions and the corresponding radiation fields. We also provide 
a first-of-its-kind glimpse at time-dependent, multi-group, multi-angle "supernova" calculations. 

To this end, we mapped a ID SESAME evolutionary calculation of a Woosley and Weaver 
(1995) 11 Mq progenitor core ~5 milliseconds after bounce into VULCAN/2D and continued the 
post-bounce calculation in full 2D for a further ~22 milliseconds. For this study, we employed only 
5 v e energy groups, with group centers at 2.5, 10.8, 30, 74.6, and 177.9 MeV. In addition, we used 
the equivalent of 16 x 10 angular bins and 102 x 51 spatial zones, the latter distributed as described 
in Ott et al. (2004). The calculation covered ~69000 timesteps, with the timesteps themselves 
ranging between ~0.1 to ~0.3 microseconds. 

Figures 2 and 3 are snapshots of the Y e distribution with velocity vectors at 10.2 and 22 mil- 
liseconds into the run. The initial perturbations imposed have axial symmetry and are completely 
artificial. These figures are 240 km x 240 km and 420 km x 420 km, respectively, on a side. There is 
a slight artifact along the poles due in part to the modest angular resolution of the run and to the 
generic singularity of polar coordinates. The positions of the stalled shock wave (roughly spherical) 
are clearly seen on Figures 2 and 3 by the jumps in the behavior of the corresponding velocity 
fields. Inward pointing arrows in the unshocked region indicate infall and accretion. The slight 
moire pattern seen in the vectors is a graphical/ visual artifact and is not in the numerical data. The 
red and off-white colors depict regions of higher Y e and proton fraction, while the lowest Y e regions 
are depicted in green and blue. The characteristic trough in Y e between the shock and the inner 
core is thereby clearly indicated. The most salient and interesting feature of the flow pattern is the 
pronounced vortical motion traced by the velocity vectors between 15 and 100 kilometers radius. 
Such prompt convection is expected right after core bounce (Burrows and Fryxell 1992; Burrows, 
Hayes, and Fryxell 1995). The motion is clockwise in the positive quadrant and counterclockwise 
in the negative quadrant, with a reflected counterpart in the bottom hemisphere. The extensions 
of the high-1^ regions follow the velocity vectors. A slight top/bottom asymmetry imposed in the 
initial model has translated into a slight top/bottom flow asymmetry which persists, quite natu- 
rally, during the run. In Fig. 3, weaker vortices are seen to be developing further out. Due to the 
low spatial resolution of the grid employed, it is unlikely that the flow field has converged and only 
the grossest overall convective/overturning pattern should be considered robust. 

However, our purpose here is to determine the 2D radiation fields that result from aspherical 
matter and composition fields in time-dependent hydrodynamic supernova flows. To depict these 
we first show the electron neutrino flux spectrum vectors at 10.8 MeV. These are provided at the 
same times as those chosen for Figs. 2 and 3. Note that 10.8 MeV is roughly at the peak of the 
emergent spectrum. Figures 4 and 5 are color maps of Y e at the fiducial times of 10.2 and 22 
milliseconds, but with the flux spectrum vectors at 10.8 MeV superposed. The shock wave is in 
the whitish region (high-Y e ) just exterior to the outer red band. Its positions are more easily seen 
in Figs. 2 and 3. Figure 6 depicts the 10.8-MeV flux vector distribution at 22 milliseconds, but 
superposed on the corresponding entropy color map. 
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As expected, at the high opacities in the inner core the fluxes are very small. Elsewhere, 
the flux vectors generally point outwards, and at large distances they are more or less radially- 
oriented with magnitudes that diminish as \. However, there are many intriguing exceptions to 
this pattern. In particular, as seen in Fig. 5 and Fig. 6 at 22 milliseconds and between 15 and 80 
kilometers (a region that brackets the corresponding neutrinosphere) the 10.8 MeV fluxes generally 
point away from the region (red) of high Y e /proton fraction and higher entropies. In fact, at 30-40 
kilometers (at modest to high optical depth), the flux vectors can be 45° off the radius vector. This 
behavior is a consequence of the higher emissivities due to electron capture in the Y" e /proton-rich, 
higher-entropy regime. As Fig. 6 indicates, due to the general vortical circulation pattern set up in 
this simulation, the entropy distribution roughly follows the Y e distribution. In the thick regions, 
the fluxes tend (at least in part) to track the local emissivity, while at low optical depths in the 
semi-transparent region the specific intensities from many angles "vector" add at a given point and 
result in a smoothed partially radially streaming radiation pattern. The tendency for the radiation 
field at large distances to be more smooth and spherically-distributed than the matter in the inner 
regions is a generic feature of multi-D transport. Such an effect is not easily nor accurately captured 
with simpler transport schemes (e.g., Burrows, Hayes, and Fryxell 1995; Rampp and Janka 2002; 
Fryer and Warren 2002; Buras et al. 2003). Even so, even at 90-100 kilometers and exterior to 
the 10.8 MeV neutrinosphere, the radiation field is not radial. As Fig. 5 shows, it is enhanced at 
angles nearest the inner proton-rich environments near 50-60 kilometers. 

One measure of anisotropy is the ratio of the radial component of the flux to its magnitude 
(cos# r =F r /\F\). Figure 7 plots cos# r versus polar angle at 10.8 MeV and 22 milliseconds. As Fig. 
7 shows, this ratio at 60 km can vary by more than 40%, while at 96 km it varies by ~5%. The 
angular scale of the variation follows the angular scale of the matter asphericity; in this calculation 
at this physical time this is ~20°-30°. As Fig. 7 shows, the degree of asphericity, as measured 
by cos6* r , decreases with increasing radius, demonstrating the "smoothing" effect at low optical 
depths. 

The dependence of the flux vector and magnitude distributions on neutrino energy is revealed 
in Figs. 8 and 9. These figures are at 22 milliseconds and are similar to Fig. 5, but are for neutrino 
energies of 2.5 MeV and 30 MeV, respectively. The low-energy neutrinos are seen to assume a 
more radial pattern starting deeper in the core, while the fluxes at the higher energies (and those 
on the tail of the emergent spectrum) are strongly concentrated around the Y e contours and in the 
optically thick regions. The latter pattern is a consequence of the higher opacities at these higher 
energies. The lower opacities at 2.5 MeV cause this component to decouple deeper in and to assume 
a more radial flux (and energy density) pattern at a given radius. Not only are the first-moment 
(flux) and zeroth-moment (energy density) anisotropies in the interior a function of neutrino energy 
and radius, but the anisotropy of the emergent spectrum at infinity is a function of energy as well. 
This means that the spectrum at a given radius is generically a weak to modest function of angle 
and will reflect the anisotropy of the matter due to rotation and convective overturn interior to 
the shock and around the respective neutrinospheres. For instance, at 22 milliseconds and 60 km, 
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cos6* r varies by as much as ~40% at 10.8 MeV (as Fig. 7 indicates), but only by ~25% at 30 
MeV and by ~5% at 2.5 MeV. At 96 km, the corresponding numbers for 2.5, 10.8, and 30 MeV 
are 0.5%, 5%, and 15%, respectively, while at 125 km they are all less than 2%. Quite naturally, 
cos6> r goes to 1.0 at large distances. However, another measure of anisotropy is the variation with 
angle of the magnitude of the flux itself (\F\) and this does not flatten at large distances. Figure 
10 depicts for radii from 60 to 125 km the dependence with polar angle of the magnitude of the 
flux at 10.8 MeV and 22 milliseconds. In particular, at 125 km \F\ for 10.8 MeV varies by as much 
as 80% from 0° to 40° polar angle and at 250 km it varies by 10-15% from top to bottom. At very 
large radii (^infinity), \F\ at 10.8 MeV varies by ~15% from top to bottom. While not small, a 
much more anisotropic matter distribution, such as what rotation or a large initial core anisotropy 
can generate, is required to produce a larger degree of flux anisotropy at infinity. Nevertheless, 
the variations with angle and radius we identify here are important new components of supernova 
theory that previous simulations could not calculate. 

Figures 11 and 12 depict the associated maps for the energy- integrated zeroth- moment (total 
v e energy density) at 10.2 and 22 milliseconds, respectively. Superposed on these plots are the 
velocity vectors. These distributions, and not the flux maps, are rough guides to the distribution 
of neutrino heating thought to be important in the neutrino-driven mechanism of core-collapse 
supernova explosions (Bethe and Wilson 1985; Burrows, Hayes, and Fryxell 1995; Liebendorfer et 
al. 2001; Fryer and Warren 2002; Rampp and Janka 2002). However, a more interesting plot is Fig. 
13, which shows the relationship at 10.8 MeV between the flux vector field and the corresponding 
energy-density pattern. In this figure, the scale is 200 km x 200 km. Combined in this way we 
can see more directly the correlation between a flux and an energy-density pattern. In addition, 
though not obvious in Figs. 4 through 13, the ray effects inherent in the S n method do not manifest 
themselves for this Ng = 16 calculation until radii of ~350 kilometers, significantly exterior to the 
radius of the stalled shock. 



6. Conclusions 

We have constructed a time-dependent, multi-group, multi-angle radiation code in two spatial 
dimensions and coupled it to a hydrodynamics code. We have used the product, VULCAN/2D, 
to simulate the early post-bounce evolution of a supernova core for 22 milliseconds. The purpose 
of this exercise was to explore for the first time the mapping between the aspherical, anisotropic 
matter distributions that obtain in the convective protoneutron star and the corresponding radia- 
tion field. Highlighted were the electron neutrino flux and energy-density distributions and their 
angle dependence. We found that these radiation distributions and their anisotropics are system- 
atic functions of neutrino energy and radius and that, though the radiation field depends on the 
distributions of the physical parameters, the non-local character of the solution to the Boltzmann 
equation results in smoother radiation patterns in optically-thin regions. In addition, we found 
that the emergent spectrum should be angle-dependent. Even though the code was constructed to 
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investigate core-collapse supernovae, it can in principle be used for other astrophysical problems. 
These might include gamma-ray bursts, accretion disk evolution, and star formation. 

This paper has been a preliminary exploration and test of the new VULCAN/2D tool. More 
details on the solution method and numerical approach can be found in Livne et al. (2004). How- 
ever, this study is but the beginning of a program of development and calculation which will include 
the incorporation of further velocity-dependent terms, the upgrade of the code to handle energy 
redistribution in an explicit manner, full simulations with the MGFLD variant, simulations with 
rotation (cf. Ott et al. 2004), and simulations with more energy groups and neutrino species. Our 
philosophy has been to develop in systematic fashion, and with a realistic timetable, an evolvingly 
sophisticated computational capability in multi-dimensional core-collapse supernova studies, in par- 
ticular, and in astrophysical radiation hydrodynamics, in general. Importantly, we are gearing up to 
address the neutrino-driven mechanism of core-collapse supernovae by calculating for much longer 
times, for different progenitors, and with more energy groups and neutrino species. The code's 
efficient parallelization in group and species implies that higher-energy-resolution simulations incur 
little penalty. Finally, since the code follows transport in the angular, as well as the radial direc- 
tions, it can explore in numerical detail and for the first time the viabilility of neutron- finger or 
doubly-diffusive instabilities in the core-collapse context (Mayle and Wilson 1988). 
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Fig. 1. — Coordinate system used for the axisymmetric transport. The radiation direction vector fi 
is defined in terms of 9 and <f>, where 9 is the angle with respect to the z axis at all spatial positions 
z and r. For example, = 0, /x = 1 is along the z axis and <fi = is in the z — r plane. 

Fig. 2. — The map of the Y e distribution at 10.2 milliseconds into the immediate post-bounce 
evolution of the core of an 11 M Q progenitor from Woosley and Weaver (1995). Superposed are 
velocity vectors. The vortical motions in this low-resolution simulation are clearly delineated. The 
scale is 240 kmx240 km on a side. Green and blue depict low Y e s from ~0.1 to ~0.25 and red and 
white depict higher Y e s. The position of the stalled shock wave is indicated by the discontinuity of 
the vector field. The slight moire pattern seen in the vectors is a graphical/ visual artifact and is 
not in the numerical data. See text for details. 

Fig. 3. — Same as Fig. 2, but at 22 milliseconds. The scale is 420 kmx420 km on a side. See text 
for details. 

Fig. 4. — At 10.2 milliseconds into the simulation, the Y e distribution using the same color map as 
found in Fig. 2. The scale is 240 kmx240 km on a side. The vectors denote the electron neutrino 
flux at 10.8 MeV (in ergs cm -2 s -1 MeV" 1 , scaled for vector rendering). 

Fig. 5. — Same as Fig. 4 (Y e distribution), but at 22 milliseconds and using a scale of 420 km x 420 
km on a side. The vectors are flux vectors at 10.8 MeV. See text for a discussion of the saliant 
features of this plot. 

Fig. 6. — An entropy map of the flow at 22 milliseconds. The scale is 420 km x 420 km on a 
side. Yellow denotes low entropy regions characteristic of unshocked material, purple denotes 
intermediate entropy regions, and red denotes highest entropy regions. The vectors are the fluxes 
at 10.8 MeV, as in Fig. 5. 

Fig. 7. — The ratio (cos 9 r ) of the radial component of the flux at 10.8 MeV and 22 milliseconds to 
the corresponding magnitude of the flux versus spatial polar angle from ir/2 to —tt/2. This ratio 
is shown for various radial distances and is one measure of the anisotropy of the radiation field. 

Fig. 8. — Same as Fig. 5 at 22 milliseconds, but for 2.5 MeV electron-type neutrinos. The vectors 
trace the corresponding flux. The scale is 420 kmx420 km on a side. See text for details. 

Fig. 9. — Same as Fig. 5 at 22 milliseconds, but for 30 MeV electron-type neutrinos. The scale is 
420 km x 420 km on a side. See text for details. 

Fig. 10. — The variation with polar angle (9) of the logarithm base 10 of the magnitude of the 
neutrino flux (\F\) at 10.8 MeV at 22 milliseconds into the calculation. The polar angle varies from 
7r/2 to — 7r/2 and the flux is in units of ergs cm" 2 s" 1 MeV -1 . This quantity is shown for various 
radial distances and is another measure of the anisotropy of the radiation field. 

Fig. 11. — The total energy-integrated energy density (in ergs cm" 3 ) of the electron- type neutrinos 
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at 10.2 milliseconds into the simulation. The scale is 240 km x 240 km on a side. Superposed are 
the corresponding velocity vectors. White and red denote high energy densities, and purple, blue, 
and green denote progressively lower energy densities. Brown denotes the lowest energy densities. 

Fig. 12. — Same as Fig. 11, but at 22 milliseconds and with a scale of 420 kmx420 km on a side. 
See text for a brief discussion. 

Fig. 13. — The energy density (in ergs cm -3 MeV -1 ) at 10.8 MeV and 22 milliseconds. The scale is 
200 km x 200 km. Pink/White denotes high energy density and brown denotes low energy density. 
The sequence white, red, purple, blue, green, yellow, and brown is a sequence of decreasing order 
of magnitude. 
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